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Abstract 

The development of efficient and accurate image reconstruction algorithms is one of the 
cornerstones of computed tomography. Existing algorithms for quantitative photoacoustic to¬ 
mography currently operate in a two-stage procedure: First an inverse source problem for 
the acoustic wave propagation is solved, whereas in a second step the optical parameters are 
estimated from the result of the first step. Such an approach has several drawbacks. In 
this paper we therefore propose the use of single-stage reconstruction algorithms for quan¬ 
titative photoacoustic tomography, where the optical parameters are directly reconstructed 
from the observed acoustical data. In that context we formulate the image reconstruction 
problem of quantitative photoacoustic tomography as a single nonlinear inverse problem by 
coupling the radiative transfer equation with the acoustic wave equation. The inverse prob¬ 
lem is approached by Tikhonov regularization with a convex penalty in combination with the 
proximal gradient iteration for minimizing the Tikhonov functional. We present numerical 
results, where the proposed single-stage algorithm shows an improved reconstruction quality 
at a similar computational cost. 

Keywords. Quantitative photoacoustic tomography, stationary radiative transfer equation, 
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1 Introduction 

Photoacoustic tomography (PAT) is a recently developed medical imaging paradigm that combines 
the high spatial resolution of ultrasound imaging with the high contrast of optical imaging [7, 35, 
53, 54, 55]. Suppose a semitransparent sample is illuminated with a short pulse of electromagnetic 
energy near the visible range. Then parts of the optical energy will be absorbed inside the sample 
which causes a rapid, non-uniform increase of temperature. The increase of temperature yields 
a spatially varying thermoelastic expansion which in turn induces an acoustic pressure wave (see 
Figure 1.1). The induced acoustic pressure wave is measured outside of the object of interest, and 
mathematical algorithms are used to recover an image of the interior. 

Original (and also a lot of recent) work in PAT has been concentrated on the problem of 
reconstructing the initial pressure distribution, which has been considered as final image (see, for 
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Figure 1.1: Basic principle of PAT. A semitransparent sample is illuminated with a short 
optical pulse. Due to optical absorption and subsequent thermal expansion within the sample an 
acoustic pressure wave is induced. The acoustic pressure wave is measured outside of the sample 
and used to reconstruct an image of the interior. 


example, [1, 9, 24, 25, 26, 27, 30, 31, 37, 34, 38, 45, 51, 55]). However, the recovered pressure 
distribution only provides indirect information about the investigated object. This is due to the 
fact, that the initial pressure distribution is the product of the optical absorption coefficient and 
the spatially varying optical intensity which again indirectly depends on the tissue parameters. As 
a consequence, the initial pressure distribution only provides qualitative information about the 
tissue-relevant parameters. Quantitative photoacoustic tomography (qPAT) addresses exactly 
this issue and aims at quantitatively estimating the tissue parameters by supplementing the 
wave inversion with an inverse problem for the light propagation in tissue (see, for example, 
[2, 5, 6, 10, 15, 16, 14, 19, 36, 39, 41, 44, 46, 47, 52, 57]). 

To the best of our knowledge, apart from the very recent work [50], all existing reconstruc¬ 
tion algorithms for qPAT are currently performed via the following two-stage procedure: First, 
the measured pressure values are used to recover the initial pressure distribution caused by the 
thermal heating. In a second step, based on an appropriate light propagation model, the spa¬ 
tially varying tissue parameters are estimated from the initial pressure distribution recovered in 
the first step. However, any algorithm for solving an inverse problem requires prior knowledge 
about the parameters to be recovered as well as partial knowledge about the noise. If one solves 
qPAT via a two-stage approach, appropriate prior information for the acoustic inverse problem 
is difficult to model, because the initial pressure depends on parameters not yet recovered. This 
is particularly relevant for the case that the acoustic data can only be measured on parts of the 
boundary (limited-angle scenario), in which case the acoustic inverse problems is known to be 
severely ill-posed. Further, using a two-stage approach, only limited information about the noise 
for the optical problem is available. 

In view of such shortcomings of the two-stage approach, in this paper we propose to recover 
the optical parameters directly from the measured acoustical data via a single-stage procedure. 
We work with the stationary radiative transfer equation (RTE) as model for light propagation. 
Our simulations show improved reconstruction quality of the proposed single-stage algorithm at 
a computational cost similar to the one of existing two-stage algorithms. Obviously our single- 
stage strategy can alternatively be combined with the diffusion approximation, which has also 
frequently been used in qPAT. In the present work we use the stationary RTE since it is the more 
realistic model for light propagation in tissue. In combination with the two-stage approach, the 
RTE has previously been used for qPAT, for example, in [5, 19, 52, 47, 39, 57]. 
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1.1 Mathematical modeling of qPAT 

Throughout this paper, let Q C M. d denote a convex bounded domain with Lipschitz-boundary 
dQ, where d £ {2,3} denotes the spatial dimension. We model the optical radiation by a function 
$: SI x —» R, where <h (x, 9) is the density of photons at location x £ Q propagating in 
direction 9 £ S^” 1 . The photon density is supposed to satisfy the RTE, which reads 

9 • V X T (x, 9) + (cr (x) + p (x)) (x, 9) 

= a(x) [ k(9,9 , )^(x,9 , )d9' + q(x,9) for (x, 9) £ fl x S ^ 1 . (1.1) 

J S ' 1 - 1 

Here a (x) is the scattering coefficient, p (x) is the absorption coefficient, and q (x, 9) is the photon 
source density. The scattering kernel k ( 9 , 9') describes the redistribution of velocity directions of 
scattered photons due to interaction with the background. The stationary RTE (1.1) is commonly 
considered as a very accurate model for light transport in tissue (see, for example, [3, 1 , 1, 33]). 

In order to obtain a well-posed problem one has to impose appropriate boundary conditions. 
For that purpose it is convenient to split the boundary T := <912 x into inflow and outflow 
boundaries, 

T_ := j(x,0) €dfixS d_1 : v(x) ■ 9 < o} , 

T + := j(x, 9) £ <912 x § d_1 : v(x) ■ 9 > 0 j , 

with z/(x) denoting the outward pointing unit normal at x £ 912. We then augment (1.1) by the 
inflow boundary conditions 

<b|r_ = / for some /: T_ —> R. (1.2) 

Under physically reasonable assumptions it can be shown that the stationary RTE (1.1) together 
with the inflow boundary conditions (1.2) is a well-posed problem. In Section 2.1 we apply a 
recent result of [20] that guarantees the well-posedness of (1.1), (1.2) even in the presence of voids 
(parts of the domain under consideration, where qi and a vanish). 

The absorption of photons causes a non-uniform heating of the tissue proportional to the total 
amount of absorbed photons, 

h{x) := n{x) / <h(x,6 ) )d6 ) for x £ H . 

Js d ~ 1 

The heating in turn induces an acoustic pressure wave p: R rf x (0,oo) —> R. The initial pressure 
distribution is given by p( ■, 0) = jh, where 7 is the Griineissen parameter describing the efficiency 
of conversion of heat into acoustic pressure. For the sake of simplicity we consider the Griineissen 
parameter to be constant, known and rescaled to one. We further assume the speed of sound 
to be constant and also rescaled to one. The photoacoustic pressure then satisfies the following 
initial value problem for the standard wave equation, 

' ( dfp(x, t) — A p(x, t) = 0 , for (x, t) £ R rf x (0,00) 

< p (x, 0) = h(x ), for x £ R d (1.3) 

„ dtp (x, 0) = 0 , for x € R d . 

The goal of qPAT is to reconstruct the parameters p and a from measurements of the acoustic 
pressure p outside H. Pressure measurements are usually taken as a function of time on parts of 
the boundary <9fl. 
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1.2 The inverse problem of qPAT 

In the following we assume that acoustic measurements are available for multiple optical source 
distributions (illuminations). For that purpose, let (<&, _/)) for i = 1,... N be given pairs of source 
patterns and boundary light sources. We use T, to denote the operator that takes the pair (/r, a) 
to the solution of the stationary RTE (1.1), (1-2) with qi and /,; in place of q and /, and denote 
by 

Hj(/x, a) (x) := /r(x) / T.;(^, er)(x, 9)d9 for x € P 

J S '*- 1 

the operator describing the corresponding thermal heating. Further, we write Wo ,a for the 
operator that maps the initial data h to the solution Wq,a h := p|,9Qx(0,oo) °f the wave equation 
(1.3) restricted to the boundary dd. Appropriate functional analytic frameworks for Tj, H, and 
Wo ,a will be given in Section 2, where we also study properties of these mapping. 

The reconstruction problem of qPAT with multiple illuminations can be written in the form 
of a nonlinear inverse problem, 


Vi = (Wo, A o Hi) (//, a*) + Zi for i = 1,..., N . (1.4) 

Here v t are the measured noisy data, the operators Wo, a ° Hj model the forward problem of 
qPAT, Zi are the noise in the data, and /r*,<7* are the true parameters. The aim is to estimate 
the parameter pair (/r*,<7*) from given data Vi, and hence solving the inverse problem (1.4). 

1.3 Outline of the paper 

In this paper we address the inverse problem (1.4) by Tikhonov regularization, 

1 N 

X E II("Wo,A o Hi) (n,a) - Vi II 2 + xn(fi,a) min, 

where 1Z is a convex penalty and A > 0 is the regularization parameter. We show that Tikhonov 
regularization applied to single-stage qPAT is well-posed and convergent; see Theorem 3.2. For 
that purpose we derive regularity results for the heating operators Hi in Section 2. To establish 
such properties we use results for the stationary RTE derived recently in [ ]. 

For numerically minimizing the Tikhonov functional we apply the proximal gradient algorithm 
(also named forward backward splitting); see Section 3.4. The proximal gradient algorithm, 
is an iterative scheme for minimizing functionals that can be written as the sum of a smooth 
and a convex part [13, 12]. For the classical two-stage approach in qPAT in combination with 
the diffusion approximation, the proximal gradient algorithm has recently been applied in [59]. 
Numerical results using the proximal gradient algorithm applied to our single-stage approach 
are presented in Section 4, where we also include a comparison with the two-stage approach. 
Of course, our single-stage approach can be combined with classical gradient or Newton-type 
schemes. The proximal gradient algorithm is our method of choice, since its is very flexible and 
fast, and can be applied for a large class of smooth or non-smooth penalties. 

2 Analysis of the direct problem of qPAT 

Before actually studying the inverse problem of qPAT we first make sure that the forward problem 
is well-posed in suitable spaces and that the data depend continuously on the parameters we 
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intend to reconstruct. For that purpose we review a recent existence and uniqueness result for 
the stationary RTE allowing for voids in the domain of interest [20]. The use of a-priori estimates 
will lead to differentiability results for the operators T,: and Hj. 

2.1 The stationary RTE 

The stationary RTE has been studied in various contexts. The most prominent, apart from the 
transport of radiation in a scattering media, is reactor physics, where the equation is used in 
the group velocity approximation of the neutron transport problem. An extensive collection of 
results regarding applications as well as existence and uniqueness of solutions can be found in 
[ ]. The analysis of the RTE becomes considerably more involved if internal voids, i.e. regions 

where scattering and absorption coefficient become zero, are allowed. 

Suppose 1 < p < oo, and denote by L P (T_, \u ■ 9\) the space of all measurable functions / 
defined on T_ for which 


ll/II.Li>(r_,|i/-0|) :— 



if p < oo 
if p = oo 


is finite. We write W p (ft x 8> d x ) for the space of all measurable functions defined on Q x 1 
such that 


$ 


up 

HWPffixS'*- 1 ) 


<f> 


IIP 

iiLr(nxs d - 1 ) 


+ ||0 -V x T 


up 

HLP(nxS d - 1 ) 


+ <b|r 


iip 

llLr(r_,M|) 


is well defined and finite (with the usual modification for p = oo). The subspace of all <F € 
W p (yi x § rf_1 ) with <h| r _ = 0 will be denoted by fU p (if x § rf_1 ). Further, for a given scattering 
kernel k £ L°°(§ rf_1 x § rf_1 ) we write K: L p (hl x § rf_1 ) —> L p (£l x § rf_1 ) for the corresponding 
scattering operator, 


(K &)(x,6)= [ k{9,9')$(x,9')d9' for (x, 9) € 0 x S'*" 1 . 

J S ^- 1 


Throughout this article, the scattering kernel k is supposed to be symmetric and nonnegative, and 
to satisfy fgd-i k (9,9') d 9' = 1 for all 9 € § rf_1 . This reflects the fact that k (■, 9') is a probability 
distribution describing the redistributions of velocity directions due to interaction of the photons 
with the background. Under these assumption, the scattering operator K is easily seen to be 
linear and bounded. 

Using the notation just introduced, the stationary RTE (1.1), (1.2) can be written in the 
compact form 

f (<9 • Vs + {n + a - crK)) $ = q in 0 x S d_1 

\ ^|r_ = / on T_ . 

By definition, a solution of the stationary RTE (1.1), (1.2) in W p is any function <F £ W p (£lx S d_1 ) 
satisfying (2.1). The following theorem, which has been derived very recently in [22], states 
that under physically reasonable assumptions there exists exactly one such solution, that further 
continuously depends on the source pattern and the boundary light source. 


Theorem 2.1 (Existence and uniqueness of solutions in W p ). LetjL,a denote positive constants, 
let fi, a be measurable functions satisfying 0 < p < Ji and 0 < a < a, and let 1 < p < oo. Then, 
for any source pattern q € L p {Vt x S^ 1 ) and any boundary light source f € L p (T_, \v ■ 9\), the 
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stationary RTE (1.1), (1.2) admits a unique solution $ € TT P (12 x S> d x ). Moreover, there exists 
a constant C p (Jl,W) only depending on p, ~p and a, such that the following a-priori estimate holds 



Proof. See [22], □ 

2.2 The parameter-to-solution operator T for the stationary RTE 

Throughout this subsection, let 1 < p < oo, and let q € L°°(12 x § d_1 ) and / € L°° (T_, \ v • 9\) be 
given source pattern and boundary light source, respectively. Further, for fixed positive numbers 
Jl,a>0we denote 

V p := |€ L p (11) x L p (12 x S^” 1 ) : 0 < a < a and 0 < fi < /izj . (2.3) 

Then T> p is a closed, bounded and convex subset of LP (12) x L p (12 x , that has empty interior 
in the case that p < oo. 

Definition 2.2 (Parameter-to-solution operator for the stationary RTE). The parameter-to- 
solution operator for the stationary RTE is defined by 

T: V p -)■ W p (n x S^ 1 ): (fi, a) ^ $ , (2.4) 

where denotes the unique solution of (1.1), (1.2). 

According to Theorem 2.1 the operator T is well defined. Note further, that T depends on 
p, q, f, Jl and a. Since these parameters will be fixed in the following and in order to keep the 
notation simple we will not indicate the dependence of T on these parameter explicitly. 

Now we are in the position to state continuity properties of T derived in [ ]. We include 

a short proof of these results as its understanding is very useful for the derivation of similar 
properties of the operator describing the heating that we investigate in the following subsection. 

Theorem 2.3 (Lipschitz continuity and weak continuity of T). 

(a) The operator T is Lipschitz-continuous. 

(b) If 1 < p < oo, then T is sequentially weakly continuous. 

Proof, (a) Let (fi, a), ( fi , a) € V p be two given pairs of absorption and scattering coefficients and 
denote by <f> := T(^,cr) and $ := T(/2,<r) the corresponding solutions of the stationary RTE. 
Since q € L°°( 12 x § d_1 ) and / € L°°(T_, \u ■ 6 1), Theorem 2.1 implies that the difference <E> — 
is an element of Wff °(12 x Further, this difference is easily seen to satisfy 

(0 ■ V* T fi + a - <tK) (6 - $) = (fi - fi) $ + (a - a)& -(<r-a)K$. 

Because K is a bounded linear operator on L p (12 x S^ 1 ), the right hand side in the above equation 
is actually contained in L p ( 12 x § d_1 ). Therefore, a further application of Theorem 2.1 yields 

11^ ~ ^*11 wuynxS^ -1 ) ^ C p (n, (T )||^||L o °(nx§ d - 1 ) ( II T ~ All lp(Q) + III _ K|| p ||cr — d|| LP (Q xS d-i)^ , 

where I denotes the identity and || ■ || p the operator norm on L p (12 x S^ 1 ). Since || < h||z / oo(Q X gd-i) 
is bounded independently of 4>, this implies the Lipschitz continuity of T. 
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(b) Let (p n , cr n )neN be a sequence in V p that converges weakly to the pair (p,a) E V p , and 
denote by <h n = T (p n ,a n ) and 4> = T(/x, cx) the corresponding solutions of the stationary RTE. 
As in (a), one argues that the difference ‘Ln — <3? is contained in x S^ 1 ) and satisfies 

(6 ■ W x + p + a - crK) ($„ -$) = (//- p n ) $ n + (a - a n )$n - (a - a n ) K<L n . 


Now, from Theorem 2.1 it follows that 6 ■ V x + p + a — a K is invertible as an operator from 
Wq(£1 x S d_1 ) to L P (Q x § d_1 ). Consequently, the inverse mapping (0 ■ V x + p + a — crK ) -1 
is linear and bounded and in particular weakly continuous. It therefore remains to show that 
(p — p n ) & n + (cr — cr n )<3? n — (a — a n ) K<f> n weakly converges to zero in L P (Q x § d_1 ). To see this, 
denote by p* = p/(p— 1) the dual index and let ip E L P *(Q x § d_1 ) be any element in the dual of 
L p (0 x S^ 1 ). By Fubini’s theorem we have 

/ (p(x) - p n (x)) $ n (x,0)ip(x,6) d(x,0) = / (p(x) - p n (x)) ( / $ n (x, 6)tp(x, 6) dd) dx . 

Jnx s d -! Jn \J s^- 1 / 


The averaging lemma (see, for example, [40]) implies that the averaging operator A: W P *{VL x 
S d_1 ) —>■ L P *(Q): <f> i->- f S d-i $(-,0)d0 is compact for 1 < p* < oo. Since (<h n ) ne N is bounded 
in W°°{Fl x § d_1 ) C W P *{VL x § d_1 ), this implies that fgd-i 4> n ( ■, 0)<p( ■, 0)dO converges to 
f S d-i 4>( •, 0)<p( •, 0)d9 with respect to || • || As p n —p we can conclude that (p — p n )$ n 

converges to zero weakly. In the same manner one shows (a — a n ) <3? n —*• 0. Finally, the equality 


/ 

Jn 


nxs d ~ 


(p - p n ) (x)(K$ n )(x,0)<p(x,0)d(x,0) = / {p-p n )(x) 


in 


<h n (x, 0)(Kip)(x, 6)d6dx 


and the use of similar arguments show that (a — a n ) K<h n —*• 0. 


□ 


For the solution of the inverse problem of qPAT we will make use the derivative of T that 
we compute next. For that purpose we call h € L p (Fl) x L p {fl x § rf-1 ) a feasible direction at 
(p, a) E V p if there exists some e > 0 such that (p, a) + eh E V p . Due to the convexity of V p we 
have (p, a) + sh E V p for all 0 < s < e. The set of all feasible directions at (p, a) will be denoted 
by V p (p,a). One immediately sees that 

V p (p, a) = L P (Q) x L P (Q x S^ 1 ) if 0 < p < p and 0 < a <a. 


For (p,cr) E V p and any feasible direction h E T> p (p,cr) we denote the one-sided directional 
derivative of T at (p, a) in direction h by 


T'(p,a){h) := lim 


T((p,a) + sh) 
s 


T{p,a) 


(2.5) 


provided that the limit on the right hand side of (2.5) exists. If both limits T'(p,a)(h) and 
T '(p,a)(—h) exist and h T'(p,a)(h) is bounded and linear, we say that T is Gataux differen¬ 
tiable at ( p,a) and call T'(p,a) the Gataux derivative of T at ( p,a ). 

Theorem 2.4 (Differentiability of T). For any ( p,a ) E V p , the one-sided directional derivative 
of T at ( p,a ) in direction (h p ,h a ) E V p (p,a) exists. Further, we have T'(p,a)(h p ,h a ) = T, 
where T is the unique solution of 

U0-V x + (p + a-aK))y = -(h lx + h a -h a K)T(p,cj) in Q x S ^ 1 

\ ’L|r_ =0 on F. 

If 0 < p <fl and 0 < a < a, then T is Gateaux differentiable at (p,c r). 
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Proof. Suppose (/q cr) £ V p and let h = £ D p (/qcr) be any feasible direction. For 

sufficiently small s > 0 write <f> s := T((/q cr) + sh) and $ := T(/q cr). As in the proof of Theorem 
2.3 one shows that := (<f> s — <F)/s is contained in Wq(U x S^ 1 ) and solves the equation 
(6- V^i/i+CT-ffK)$ s = — {h p +h a — h a K)$ s . Consequently the difference — T € lTg(ffx§ d_1 ) 
solves 

(0 • V* + q + u - uK) (T s - T) = -(h p + h a - /i ct K)($ s - $). 

Application of the a-priori estimate of Theorem 2.1 shows the inequality ||T S — v f'||iyp(nxS d - 1 ) — 
Cp(%v)\\®s - ^’lli,°°(f2xS d - 1 )(ll / VlUp(n) + H^o-llLr^xS^- 1 )) Together with the continuity of T this 
implies that the one-sided directional derivative T '(/i, cr)(h ) exists and is given by lim s _^o = T. 
Finally, if 0 < q < ft and 0 < a < cf, then h i->- T'(q> a)(h) is bounded and linear and therefore T 
is Gateaux differentiable at (q, cr). □ 

Note that for any parameter pair (q, a) £ V p , the solution of (2.6) depends linearly and 
continuously on (/qV<r) £ L p { fl) x L p (Pl x § d_1 ). As a consequence, the one-sided directional 
derivative can be extended to a bounded linear operator 

TV, a): L p (Ll) x L p (Ll x S^ 1 ) -> W p (tt x S'* -1 ): (h p , h a ) ^ T , (2.7) 

where T is the unique solution of (2.6). We refer to this extension as the derivative of T at (q, a). 

2.3 The operator H describing the heating 

Throughout this subsection, let q £ L°°(Q x S d_1 ) and / £ L°° (T_, \ u ■ 0|) be given source pattern 
and boundary light source, respectively. As already mentioned in the introduction, photoacoustic 
signal generation due to the absorption of light is described by the operator 

H: V p ^L p (Q): (q,c t) ^ q [ T(q, <r)( •, 6)d9 . 

JSd - 1 

To shorten the notation, in the following we will make use of the averaging operator A: W p (Ii x 
gd-i) £P(fJ) defined by A<F = f S d-i < f ) ( ■, 0)d9. By Holders inequality the averaging operator is 
well defined, linear and bounded. Using the averaging operator we can write H(q, a) = qAT(q, a). 

Theorem 2.5 (Lipschitz continuity and weak continuity of H). 

(a) The operator H is Lipschitz continuous. 

(b) If 1 < p < oo, then H is sequentially weakly continuous. 

Proof, (a) Suppose that (q,cr), (A,d) € T)p are two pairs of admissible absorption and scattering 
coefficients. The decomposition H(q, a) = /iAT(q, cr) and the triangle inequality imply 

IIA*AT(/x, a) - AAT(A,<j)|| LP(n) 

= HMTQq a) - AAT(/i, cr) + AAT(q, a) - AAT(A, <t)|| iP(n) 

— II AT(g, cr) || L \\h ~ All lp(Q.) T || All l°°(Q) II AT(/i, cr) — AT (A, (?) || • 

According to Theorem 2.3, the operator T is Lipschitz continuous. Because A is linear and 
bounded, also the composition AT is Lipschitz. Noting that ||AT(q,cr)||l/ x, (n) and HAlIl, 00 ^) are 
bounded by constants independent of q, cr and ft, <7, this implies the Lipschitz continuity of H. 


(b) Let (/i n ,<7 n )ne n be a sequence in V p that converges weakly to (fa, a) € V p . Since T is 
weakly continuous and A is linear and bounded, (AT(/j n ,ff n )) n6 pj converges weakly to AT(fi, a). 
Further, for any function (p € L P *(Q), the dual space of L P (D), we have 


(n(x)(AT)(fa,a)(x) - /i„(i)(AT )(fi n ,a n )(x)) ip(x) dx 


— II AT(/x, cr) ||x<x>(n) 
+ II Ahi II 


in 


(h(x) - /j, n (x)) ip(x) dx 
((AT)(/x, cr)(x) - (AT)(// n , cr n )(x)) (p(x)dx 


The weak convergence of fi n and (AT(/i„,ff n )) ne N therefore implies the weak convergence of 
H n AT(/j, n ,a n ) to /tAT(/t, cr) and shows the weak continuity of H. □ 


Note that for the case 1 < p < oo, the averaging operator A is even compact (see [40]) which 
implies the compactness of the composition AT. As a consequence, for any given /t, the partial 
mapping a i->- g(AT)(/i, a) is compact. It seems unlikely, however, that the full operator H is 
compact, too. 

Theorem 2.6 (Differentiability of H). For any (fi,cr) € T> p , the one-sided directional derivative 
of H at (fi,cr) in any feasible direction (h p ,h a ) € V p (fi,cr) exists. Further, we have 


H '(fi,a)(h tl ,h tT ) = h p [ T(fi,a)(-,Q)d9 + fi f T'(fi,a)(h p ,h a )(-,9)d9, (2.8) 

JS d ~ 1 JS d ~ 1 

where T'(/t, cr)(h p , h a ) denotes the one-sides directional derivative of T at (fa, a) in direction 
(/t M ,/i CT ) and can be computed as the solution of (2.6). Finally, ifO<p<Ji and 0 < a <W, then 
H is Gateaux differentiable at (fa, a). 

Proof. Let (fa, cr) € V p and let (h p ,h a ) € V p (fa,a) be a feasible direction. For sufficiently small 
s > 0, we have 


H((/t, cr) + s(h p , h a )) - H(/i, a) 


(fa + sh p )( AT) ((fa, a) + s(h p ,h a )) - //(AT) (fa, a) 


= hfj( AT) ((fa, a) + s(h p , h a )) + fa 


(AT) ((//,<cr) + s(VM) - ( AT ) (h’V) 


According to Theorem 2.3, the operator T is continuous and therefore the first term converges to 
h p (AT)(/t, a) as s —>• 0. Because T is one-sided differentiable, see Theorem 2.4, the second term 
converges to //AT'(g, cr)(h p , h a ). Finally, if 0 < fa < Ji and 0 < a < a, then H'(fa, a)(h) is linear 
and bounded in the argument h which implies the Gateaux differentiability of H at (fa, a). □ 


Recall that for any (n, a) € V p , the derivative T '(fa, a) is bounded and linear. Therefore, the 
right hand side of (2.8) depends linearly and continuously on (h p , h a ) € L p ( fl) x L P (Q. x S^ -1 ). As 
a consequence, the one-sided directional derivative of H at (/t, cr) can be extended to a bounded 
linear operator H'(/i, <r): L p (0) x L P (Q. x § rf_1 ) —>• L P (D). We will refer to this extension as 
the derivative of H at (fa, a). The derivative of H can be written in the form H '(fa,a)(h) = 
h p A(T(/i, cr)) + fa A (T'(fa, a)(h)). 
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2.4 The wave operator W^,a 

Let U C M. d be a bounded and convex domain with smooth boundary. We assume that Q C U 
and write Lf l (R d ) for the space of all square integrable functions defined on M. d that are supported 
in Q. Likewise we denote by the space of all infinitely differentiable functions defined on 

having support in fL Further, let A C dU be a relatively open subset of dU, and denote by 
diam({7) the maximal diameter of U and by dist(f2, A) the distance between and the observation 
surface A. 

Definition 2.7 (The wave operator Wjj a)- Let ico.a : (0, oo) -> 1 be a smooth, nonnegative, 
compactly supported function with W(i,A(t) = 1 for all dist(fl, A) < t < diam(Z7) — dist(0, A). We 
then define the wave operator by 

Wf!, A : C(?(R d ) C Li(R d ) —)• L 2 (A x (0,oo)) : h ^ w q ,ap\ax( 0 ,oo) , (2-9) 

where p denotes the unique solution of (1.3). 

The operator W^ ]A maps the initial data of the wave equation (1.3) to its solution restricted 
to A C dU and models the acoustic part of the forward problem of PAT. The cutoff function 
A accounts for the fact, that in the two dimensional case the solution of the wave equation has 
unbounded support in time but measurements can only be made over a finite time interval. 

In the following we use a result from [42] to show that A is a bounded linear and densely 
defined operator, and therefore can be extended to a bounded linear operator on Lq(M. d ) in a 
unique manner. 

Theorem 2.8 (Continuity of the wave operator Wfi |A ). There exists some constant cn, A such that 
||Wn iA /i||i2( Ax (o,oo)) ^ Cfi iA ||/i ||£2 ( K d> for all h G C'^’(M‘ i ). Consequently, there exists a unique 
bounded linear extension 


"WY 2 iA : Li(M. d ) —» L 2 (A x (0, oo)) with W^aIc^^) — Wn jA . 


With some abuse of notation we again write Wq )A for Wq jA in the sequel. 


Proof. It is sufficient to consider the case when the data are measured on the whole boundary 
A = dU. The well known explicit formulas for the solution of (1.3) in two and three spatial 
dimensions (see, for example, [23, 32]) imply that for every (y,t) G dU x (0, oo) we have 


(W n , A h) ( : y,t ) 


9t £ f S d- 1 h(y + ru) dojdr for d = 2 

d t (t / sd _! h(y + tu) dw) for d = 3 . 


( 2 . 10 ) 


We define the spherical mean Radon transform M: Cq^M^) —»• C°°(dU x (0, oo)) by ML (y, t ) := 
1/ud-i fgd-i h(y + tuj)duj for (y,t) € dU x (0, oo). With the spherical mean Radon transform, the 
wave operator can be written as (Wn ,\h) (y, t ) = wn,A(t) dt Jq rM/i (y, r) /Vt 2 — r 2 dr in the case 
of two spatial dimensions and (Wg A h) ( y,t ) = wn,A(t) dt (fM/j) (y, t) in the three dimensional 
case. 

Next we use a Sobolev estimate derived in [42], which states that for every A G R there 
exists a constant ck,\ such that ||M/i|| ff A+(d-i)/ 2 (g r / x (o !00 )) < ck ,a ll^ll_ffA(r/) f° r an Y h G Cff(W l ). 
Application of this identity with A = 0 and using the smoothing properties of the Abel transform 
by degree 1/2 for the case of two spatial dimensions yields the continuity of W^ iA with respect 
to the L 2 topologies. In particular, W^ iA has a unique bounded linear extension to Li(M. d ). □ 
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For solving the inverse problem of qPAT we will further utilize au explicit expression for the 
adjoint of Wn,Aj that we compute next. 

Proposition 2.9 (Adjoint of the wave operator). For v £ L 2 ( A x (0, oo)) n C 1 (A x (0, oo)) ; and 
every x £ P, we have 


(W* nA v) (x) = 


J _ f f°° dt{wu,Av) ( y,t ) 

2?r Ja J\x—y\ y/r 2 - \x - y\ 2 

1 f d t (wn t Av)(y, \x-y\) 


drdS(y) if d = 2 


( 2 . 11 ) 


An 


'A 


\x-y\ 


d S(y) 


if d = 3 . 


Proof. This is a simple application of Fubini’s theorem and the explicit expression for Wn,A h 
given in (2.10). □ 


3 Single-stage approach to qPAT 

In this section we solve the inverse problem of qPAT by a single-stage approach. Our setting 
allows acoustic measurement for multiple sources. Such a strategy has been called multi-source 
qPAT or multiple illumination qPAT (see [(i, 17, 58]). For that purpose, throughout this section 
qt € L°°(P x § d_1 ) and fi € L°° (r_, \v ■ (9|), for i = 1,... IV, denote given source patterns and 
boundary light sources, respectively. Recall that P C denotes a bounded convex domain with 
Lipschitz boundary and T_ denotes the inflow boundary consisting of all pairs (x, 9) £ SP x 
with v(x) ■ 9 < 0. 


Figure 3.1: Setup for single-stage qPAT. 
The stationary RTE governs the light propaga¬ 
tion in the domain P. the absorption of photons 
induces an initial pressure wave proportional to 
the heating Hj(/x, a). Further, P is supposed to 
be contained in another domain U, and the pres¬ 
sure waves are measured with acoustic detectors 
located on a open subset A C dU of the bound¬ 
ary of U. 

To indicate the dependence of the solution of the stationary RTE on the pair ( qi , fi) we write 
Tj: T >2 —>• L 2 (P x § d_1 ) for the solution operator of the stationary RTE (1.1), (1.2) with sources 
(qi, fi). Here V 2 is the set of all admissible pairs (fi, a) defined in (2.3). Further, we use 

Hi: V 2 L 2 (P): (fi, a) ^ y [ T,(y, a)( ■, 9)d9 

to denote the corresponding operator describing the heating. For the coupling to the acoustic 
problem, it will be convenient to consider Hj(/r, a) £ L 2 (P) as an element of L^(M rf ), by extending 
it to a function defined on W l that is equal to zero on M. d \ P. 

Further, recall the Definition 2.7 of the wave operator Wjj^ modeling the acoustic problem, 
that maps the initial pressure h in the wave equation (1.3) to its solution restricted to A C dU. 
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Here U is a convex domain with smooth boundary that contains the support of /. To apply the 
results of Section 2, in the following we assume that Cl C U. Then, according to the Section 2, the 
operators H* are Lipschitz continuous, weakly continuous and one-sided directional differentiable, 
and WgA is linear and bounded. A practical representation of these domains is illustrated in 
Figure 3.1. 

3.1 Formulation as operator equation 

In order to apply standard techniques for the solution of inverse problems we write the recon¬ 
struction problem of (multiple-source) qPAT as a single operator equation. For that propose we 
denote by 

F: V 2 ^ (L 2 (A X (0, oo))) N 

(->• (Wn, A ° Hi(pt, cr),..., Wfi >A o H N (p,cr)) 

the operator describing the entire forward problem of qPAT. Further we denote by ||u||^ := 
Y^iLi 11 11 z . 2 (Ax (o,oo)) the s q uare d norm on L 2 ( A x (Cfoo))^. 

Theorem 3.1 (Properties of the forward operator of qPAT). 

(a) The operator F is sequentially weakly continuous. 

(b) The operator F is Lipschitz continuous. 

(c) For any (p,cr) € D 2 , the one-sided directional derivative in any feasible direction h € 
T >2 (yu, ct) exists. Further, 

F \n,o)(h) = (Wn,A ° Hi(//, cr)(It),..., Wn jA o tl' N (p,a)(h)) (3.1) 

where H((ju, cr)(h) is given by (2.8) with T replaced by T 

(d) If 0 < p < Ji and 0 < a < 7f, then F is Gateaux differentiable at (p,a). 

Proof. All claims follow from the corresponding properties of the operators (see Theorems 2.5 
and 2.6) and the boundedness of WgA discussed in Theorem 2.8. □ 

The inverse problem of qPAT with multiple illuminations consists in solving the nonlinear 
equation 

v = F(ji*,<r*) + z, (3.2) 

where (//*, a *) is the unknown, v = (v \,..., vjy) are the given noisy data, F: T >2 —>• L 2 (Ax (0, oo)) N 
is the forward operator, and z is the noise in the data. Our single-stage approach for qPAT consists 
in estimating the parameter pair (/P,<7*) directly from (3.2). In contrast, existing two-stage 
approaches for qPAT first construct estimates hi for the heating functions Hj(/P, a*) from data 
Vi by numerically inverting WgA, and subsequently solve (hi,... /ijv) = (Hi(/x, a), ... , H n(h, u)) 
for (p.,a). 

There are at least two common methods for tackling an inverse problem of the form (3.2): 
Tikhonov type regularization methods on the one and iterative regularization methods on the 
other hand. In the following we apply Tikhonov regularization to the inverse problem of qPAT. 

3.2 Tikhonov regularization for single-stage qPAT 

We address the inverse problem (3.2) by Tikhonov regularization with general convex penalty. For 
that purpose, let TZ: L 2 (fi) x L 2 (SI x S^ 1 ) -^MU {oo} be a convex, and lower semicontinuous 
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functional with domain V(JZ) := {(p,,a) E L 2 (H) x L 2 (fl x § d_1 ) : TZ(p,,a) < oo}. We assume 
that TZ is chosen such that V 2 n V(TZ) is non-empty. 

Tikhonov regularization with penalty 1Z consists in computing a minimizer of the generalized 
Tikhonov functional 

T v ,\: L 2 (fi) xL 2 (Hx S^" 1 ) -4lu { 00 } 

. f±\\F(fx,a)-v\\ 2 N + \K(ii,a) if (», a) € V 2 n V(K) (3.3) 

I 00 otherwise. 

Here A > 0 is the so called regularization parameter which has to be chosen accordingly, to balance 
between stability with respect to noise and accuracy in the case of exact data. The data-fidelity 
term i||F(/i, a) — v\\^ guarantees that any minimizer of (3.3) predicts the given data sufficiently 
well. The regularization term XTZ(p,, a) on the other hand avoids over-fitting of the data and 
makes the reconstruction process well-defined and stable. 

Note that Tikhonov regularization with penalty 1Z is designed to stably approximate a solution 
of the constrained optimization problem 

lZ(n,a) —>• min such that FHr, <r) = v*. (3-4) 

(fj,,o-)£T>2r\T>(TZ) 

Here v* E ran(F) is an element in the range of F and is referred to as exact data. Any solution of 
(3.4) is called 7£-minimizing solution of the equation F(/r, a) = v*. Under the given assumptions 
there exists at least one 7£-minimizing solution, which however is not necessarily unique, see 
[18, 49], 

The properties of the operator F derived above and the use of general results from regular¬ 
ization theory yield the following result. 

Theorem 3.2 (Well-posedness and convergence of Tikhonov regularization). 

(a) For data v E L 2 ( A x ( 0 , 00 ))^ and every A > 0, the Tikhonov functional T V} \ has at least 
one minimizer. 

(b) Let A > 0, v E L 2 { A x ( 0 , 00 ))^, and let (v n )neN &e a sequence in L 2 ( A x ( 0 , 00 ))^ 
with ||u — v n \\ n —>• 0. Then every sequence of minimizers (/a n ,a n ) E argminT^A has a 
weakly convergent subsequence. Further, the limit u of every weekly convergent subsequence 
(jU r (ra),°V(n))neN is a minimizer (/u,, cr) of T\ v and satisfies 7£(// r ( n ), a r(n )) ->• lZ(fi,a) for 
n —>• 00 . 

(c) Let v * E ran(F), let (<5 n ) ngN C (0,oo) be a sequence converging to zero, and let (t'n)nGN C V 
be a sequence of data with ||u* — v n \\ N < h n . Suppose further that (A n ) nS N C (0, 00 ) satisfies 
\ n —>• 0 and 6 2 /X n —> 0 as n —>• 00 . Then the following hold: 

• Every sequence {(x n ,a n ) E argmin7), fc ,A fc has a weakly converging subsequence. 

• The limit of every weakly convergent subsequence (/h-(ra)> °Y(n))raeN °f (l^n, ^n)ne N is 

an TZ-minimizing solution of F (fi, a) = v* and satisfies LZ{p, T t n \,a T i n ^) —>• 

W,a*). 

• If the TZ-minimizing solution o/F(/x, a) = v* is unique, then (/j, n ,a n ) —*■ 

Proof. Since F is sequentially weakly continuous (see Theorem 3.1) and T> 2 is closed and convex, 
this follows from general results of Tikhonov regularization with convex penalties, see for example, 
[48, Th ru . 3.3, Thru. 3.4, T hr u. 3.5], □ 
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3.3 Gradient of the data-fidelity term 

For numerically minimizing the Tikhonov functional we will use the gradient of the data-fidelity 
term 

1 | N 

•^(/b cr ) := - ||F(/x, a) - v||^ = - ^ ||'W^aMT,^, a) ~ Vi||J, a(Ax(0>oo)) . (3.5) 

i =1 

Recall that € L 2 (A x (0, oo)) are the given data, T, is the solution operator for the stationary 
RTE with source patterns qi and boundary light sources /), A is the averaging operator, and 
W q a is the solution operator for the wave equation. 

Let (fi, a) € V 2 be some admissible pair of parameters and let (h^,h a ) e-x <r)(/i M , h a ) 

denote the one-sided directional derivative of T at (//, < 7 ). We define the gradient VJ 7 (//,cr) of T 
at (fa, a) to be any element in L 2 {fl) x L 2 (H x S d_1 ) satisfying 

{VJ r (/r,o-),(h M ,h (7 )) i2(a)><L2(nxSti _ 1) = J' / (^,cr)(/i At ,h (7 ) for {h^h a ) eV 2 (n,a). (3.6) 

From Theorem 3.1 and the chain rule, it follows that cr)(hfj,, h a ) exists for any feasible 

direction {h^,h a ) € T> 2 (p,a). Further, in the case that p and a are strictly positive, we have 
V 2 (p,a) = L 2 [fl) x L 2 {fl x S^ -1 ), which implies that \/F{p.,a) is uniquely defined by (3.6). 

In order to compute the gradient we derive a more explicit expression for the one-sided direc¬ 
tional derivative. 

Proposition 3.3 (One-sided directional derivative of the data-fidelity term). Let (p,cr) € V 2 be 
an admissible pair of parameters and let (h^fli^) € V 2 (n,a) be a feasible direction. Then we have 

N 

•F {hi ° r )(^7i, her) E (A^W^a [W^ A gATj(g, a) - Vi ] - A(^), K) 

i— 1 

N 

+ J2{-^*i+(^iWi,K) L 2^ x§d -i ) , (3.7) 
2—1 

where 4>j := Tj(/x, cr), and <3?* is the unique solution of the adjoint equation 

(-0 ■ V* + (^ + <7 - crK)) <S>* = A>W^ a [W n , A //A ^ - vf infix (3.8) 

satisfying the zero outflow boundary condition <h*|r + = 0. 

Proof. Obviously it is sufficient to consider the case N = 1, where we write v, T, and <h* in 
place of Vi, Tj, <3?j and <h*. By (3.5) we have 

J r '(^,cr)(h M ,/i (T ) 

= {Wn, A /iAT(g, cr) - v,W n ,Ah»AT(fi,a) + W n , A /iAT'(g, (t)(/i m , V)) l2(Ax(0iOo)) 

= (W a ,AMA$ - v, Wo iA h M AT) L2(Ax(0 oo)) 

+ <Wo,AhAT - v, Wn, A /iAT'(ju, <r)( V K)) L 2 (Ax(0 oo)) 

= < A$W£, a [Wn^AS - v ], K) L2{a) 

+ (A*/xWq A [Wn iA /iAT - v] ,T / (/r,cr)(h /i ,h <T )) L2(nxSd _ 1) . (3.9) 
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Recall that 4>* is the solution of (3.8) with source term q = A*/iW^ A [WqA/xAT — v] and 
zero outflow boundary conditions <h*|r + = 0. Further, according to Theorem 2.4, the one-sided 
directional derivative of T is given by T 7 (/x, h a ) = T, where 'h is the unique solution of 

(6 ■ S/ x + {fi + o — o‘K))'F = —(hp + h a — h a K)<3? with inflow boundary conditions ’F|r_ = 0. The 
zero outflow and zero inflow boundary conditions of 4>* and T, respectively, and one integration by 
parts, show (— 0- V x 4>*, 4 , )L 2 (Q X §d-i) = (4>*, 0-S7 x '&)i, 2 (nx§ d - 1 )- Further, notice that the averaging 
operator A has adjoint A*: L 2 (fl) —>• L 2 (H x § d_1 ), (A *g)(x,v) = g(x), and that the scattering 
operator K is self-adjoint. 

Using these considerations, the second term in (3.9) can be written as 

(A*juWn >A [Wa,A/iA<f> - v], T'(/x, a)(h^ ^)> L2(nxSd _ 1) 

= ((-0 • V* + ( M + a - aK)) <T, 4-)^^ 

= (4>*, {0 ■ V x + (/x + a - crK)) L^(nx§d-i) 

= , — (. hfj, + h a — h a JZ) <3?) L 2 (QxS^- 1 ) 

= - hfj, + ha) ^(Qxgd-i) + ((K<f>)4>*, /i CT ) L 2 (nxS d-i) 

= - + + ((K<h)4>*, /r 0 -) i 2 (nxS d- 1) 

= (-A(^*), h IM ) L2(n) + (-$$* + (K<&)**, Mz^nxs*-!) • 

Together with (3.9) this yields the desired identity (3.7). □ 

Let (/r, a) G V 2 be an admissible pair of absorption and scattering coefficient. If /x, u are 
both strictly positive, then one concludes from Proposition 3.3 that the gradient of T at (/x, <r) is 
uniquely defined and given by \7F(g,cr) = (V^J-^/x, <r), V CT J r (/x, a)) with 

N 

V^(/x,a) = (A^W^a [W(j iA /xA$i - Wi ] - A(4>j<h*)) 

2=1 
N 

V^(/x,a) = ^ (-Tdh* + (K<3?i)<f>*) . 

2=1 

Here <3?* := T,(/x, <r), and is the solution of the adjoint equation (3.8) with zero outflow 
boundary condition. In the case that /x, a are not both strictly positive, the gradient is not 
uniquely defined by (3.6). However, Proposition 3.3 implies that the vector V.F(/x, <r) defined 
by (3.10), (3.11) still satisfies (3.6). We therefore take (3.10), (3.11) as gradient of T at any 
Ox, 0 ") G V- 2 . 


(3.10) 

(3.11) 


3.4 Proximal gradient algorithm for single-stage qPAT 

In order to minimize the Tikhonov functional we apply the proximal gradient (or forward backward 
splitting) algorithm, which is an iterative algorithm for minimizing functionals that can be written 
as a sum T + Q, where T is smooth and Q is convex [13, 1 ]. The proximal gradient algorithm 
computes a sequence of iterates by alternating application of explicit gradient steps for the first 
functional T and implicit proximal steps for the second functional Q. 

To apply the proximal gradient algorithm for minimizing the Tikhonov functional (3.3) we 
take F(g,a) = ^||F(/x,<r) — u||)y for the first and Q(g,a) = A7£(/x,< 7 ) for the second functional. 
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The proximal gradient algorithm then generates a sequence (fJ, n ,a n ) of iterates defined by 

(fi n+ i,a n+ i) := prox SnA7e a n ) - s^VJ 7 (//„, a n )) for n € N . (3.12) 

Here (/J,o,<Jo) €. T >2 H V(1Z) is an initial guess, s n > 0 is the step size in the n-th iteration, 
V J~ (fi, a ) = (V /i J 7 (/r, cr), V CT J r (/r, cr)) is the gradient of J 7 given by (3.10), (3.11), and 

1 2 

P rox ^A^ (A. o') := argmin - \\(p, a) - {p,, a)\\ + s n \K(n,<r) (3.13) 

(<r,/i)ei>2nz>(72.) w 

is the proximity operator corresponding to the functional s n XR.(/j,,cr). 

Remark 3.4 (Lipschitz continuity of VJ 7 ). Note that the gradient VJ 7 of the data-fidelity term is 
easily shown to he Lipschitz continuous. This either can be deduced from the explicit expressions 
(3.10), (3.11) or by usingV F{p,,a) = F^/i, <j)*(F (/j,,a) — v). In any case, the Lipschitz continuity 
of (p,,a) i->- VIF(h,g) follows from similar arguments as in the proofs of the Lipschitz-continuity 
of T and H, given in Section 2.2 and Section 2.3, respectively. 

Convergence of the proximal gradient algorithm (3.12) is well known for the case that J- is 
convex with /3-Lipschitz continuous gradient and step sizes satisfying s n E [e, 2/(5 — e] for some 
constant e > 0, see [13, 12]. These results are also valid for infinite dimensional Hilbert spaces. 
Because our forward operator F is nonlinear, the data-fidelity term J 7 is non-convex and these 
results are not directly applicable to qPAT. Note however Recently, the converge analysis of the 
proximal gradient analysis has been extended to the case of non-convex functionals; see [4, 8, 11]. 


4 Numerical implementation 


Our numerical simulations are carried out in d = 2 spatial dimensions. The stationary RTE is 
solved on a square domain H = [— 1, l] 2 . For the scattering kernel we choose the two dimensional 
version of the Henyey-Greenstein kernel, 


k(9,9') 


1 1 -g 2 

2n 1 + g 2 — 2g cos(0 • 9') 


for 9, 9' € S 1 , 


where g E (0,1) is the anisotropy factor. Before we present results of our numerical simulations 
we first outline how we numerically solve the stationary RTE in two spatial dimensions. This step 
is required for evaluating both, the forward operator F and the gradient VJ 7 of the data-fidelity 
term. 


4.1 Numerical solution of the RTE 


For the numerical solution of the stationary RTE (1.1), (1.2) we employ a finite element method. 
For that purpose one calculates the weak form of equation (1.1), (1.2) by integrating the equation 
against a test function w: H x S 1 —>• K. Integrating by parts in the transport term yields 



n J s 1 


(— 9 ■ \7 x w + fiw + aw — aKw) $ d0dx + / 4>u> (9 ■ v) dcr = / / qw (19 dx . (4.1) 

Jdfix s 1 Jn J § 1 


Here we dropped all dependencies on the variables to shorten notation and dcr denotes the usual 
surface measure on dli x S 1 . 
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The numerical scheme replaces the exact solution by a linear combination in the finite element 
space 


N h 


4> ( h \x,9 ) = (x,d) 


(4.2) 


1=1 


where any ij/ j h ' 1 (x, 9) is the product of a basis function in space and a basis function in velocity 
and the sum ranges over all possible combinations. The spatial domain is triangulated uniformly 
with mesh size h as illustrated in Figure 4.1. The velocity direction on the circle is divided into 
16 equal subintervals. We use P\ Lagrangian elements, i.e. piecewise affine functions, in the two 
dimensional spatial domain as well as for the angle. 



Figure 4.1: Spatial finite element dis¬ 
cretization. The square domain fl = [—1, l] 2 
is divided into 2N 2 triangles. To any of the 
(N + l) 2 grid points x ±,..., X(jv+i )2 a piecewise 
affine basis function is associated, that takes the 
value one at one grid point and the value zero 
on all other grid points, and is affine on every 
triangle. 


To increase stability in low scattering areas we add some artificial diffusion in the transport 
direction. This is called the streamline diffusion method, see for example [33] and the references 
therein. In the streamline diffusion method the solution is approximated in the usual way by 
(4.2). However, the test functions take the form 


N h 

w(x, 9) = Wj(ipj(x, 9) + 5(x, 9) 6 ■ V x i/jj(x, 9 )), (4.3) 

l=i 

where the additional term introduces some artificial diffusion. In our experiments, the stabilization 
parameter is taken as 5(x, 9) = 3/i/100 for fj,(x ) + a(x , 9)+ < 1 and zero otherwise. Note that the 
streamline diffusion method provides a fully consistent stabilization of the original problem. 

Making the ansatz (4.2) for the numerical solution and using test functions of the form (4.3), 
equation (4.1) yields a system of linear equations for the coefficient vector of 

the numerical solution. The entries of and can be calculated by setting = i/ji and 

w = tpj + 59 ■ V x^j ■ F° r simplicity we only consider the case q = 0 corresponding to zero sources 
of internal illumination. Then, similar to (4.1) we obtain 


o J s 1 


(56 ■ X7xipi — i/ii) 9 ■ V x i/jj d0dx + / \9 ■ v\ 



n 4s 1 


(fj, + <7 — uK) ( , if>j + 59 ■ V x iJjj) if)id9dx = J \9 • u\ . (4.4) 


The entries of the matrix now can be calculated by evaluating the integrals on the left hand 
side of (4.4). The right hand side of (4.4) together with the prescribed boundary light sources on 
T_ yields the entries of b^ h \ 
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Figure 4.2: Original phantom and simulated data. True absorption coefficient (left), simu¬ 
lated heating function (middle), and simulated pressure data on the boundary (right; the detector 
location varies in the horizontal and time in the vertical direction). 


4.2 Numerical results 

For the following numerical results, the stationary RTE is solved by the finite element method 
outlined in Subsection 4.1. For that purpose the domain = [—1,1 ] 2 is discretized by a mesh of 
triangular elements (compare Figure 4.1). The angular domain is divided into 16 subintervals of 
equal length. The anisotropy factor is taken as g = 0.6 and the scattering coefficient is taken as 
<7 = 3. We use a single boundary source distribution / representing a planar illumination along the 
lower edge [—1,1] x {—1}, where all photons enter the domain vertically. For solving the inverse 
problem we used a mesh containing 7442 triangular elements. In order to avoid performing inverse 
crime, for simulating the data we used a finer mesh containing 20402 triangular elements. 

The solution of the two-dimensional wave equation (1.3) is computed by numerically evaluat¬ 
ing the solution formula (2.10), where the detection curve A = {3/2(cos ip, sin p) : p E (— it, 0)} is 
a half-circle on the boundary of - 63 / 2 ( 0 ). The adjoint is evaluated by numerically imple¬ 

menting (2.11). This can be done efficiently by a filtered backprojection algorithm as described 
in [9, 25]. The geometry of 12 and A is similar to the one illustrated in Figure 3.1. 




0.9 



Figure 4.3: Reconstruction results for SIMULATED DATA. Reconstructed absorption coef¬ 
ficient using our single-stage approach (left), and reconstructed absorption coefficient using the 
usual two-stage approach (right). 

For our initial experiments we assume the scattering coefficient a to be known. In such a 
situation, the proximal gradient algorithm outlined in the Subsection 3.4 reads 

Pn+I ■■= prox SnA7e (p n - SnV^J 7 (//„, cr)) for n € N 
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where s n > 0 is the step size, V M .F(/r, a) is the gradient of T in the first component, given by (3.10), 
and prox SnA 7 £ (•) is the proximity operator similar as in (3.13). In the presented numerical exam¬ 
ples the regularization term is taken as a quadratic functional TZ(fi) = \ \\dx^\\i, 2 (^ ) + \ \\d y n\\ 2 L 2 ^y 
In order to speed up the iterative scheme we compute the proximity operator only approximately 
by projecting the unconstrained minimizer argmin^ ^\\fi — £l\\ 2 + s n XlZ(ij) on V 2 . Therefore the 
main numerical cost in the proximal step is the solution of a linear equation, which is relatively 
cheap compared to the evaluation of the gradient V^J 7 . 

In Figures 4.2 and 4.3 we present results of our numerical experiments. Figure 4.2 shows 
the true absorption coefficient as well as the heating function and the simulated pressure data. 
The left image in Figure 4.3 shows the numerical reconstruction with the proposed single-stage 
approach using 40 iterations of the proximal gradient algorithm. We observed empirically, that 
the proximal gradient algorithm stagnated after 20 to 40 iterations and therefore we used 40 as a 
stopping index. For comparison purpose, the right image in Figure 4.3 shows reconstruction results 
using the classical two-stage approach. For that purpose we apply Tikhonov regularization and the 
proximal gradient algorithm to the inverse problem h = Hj(^r) + z^. Again the iteration has been 
stopped after 40 iteration, where the iteration has been stagnated. The approximate heating 
h is computed numerically by applying the two dimensional universal backprojection formula 
[9, 38, 29] to the acoustic data v = W^,a ° Hj(/r) + z. All computations have been performed in 
Matlab on a notebook with 2.3 GHz Intel Core i7 processor. The total computation times have 
been 26 minutes for the two-stage approach and 38 minutes for the single-stage approach. 

One notices that in both reconstructions some boundaries in the upper half are blurred. Such 
artifacts are expected and arise from the ill-posedness of the acoustical problem when using 
limited-angle data; see [28, 43, 56]. However these artifacts are less severe for the single-stage 
algorithm than for the classical two-stage algorithm. Further, in this example, the single-stage 
algorithm also yields a better quantitative estimation of the values of the absorption coefficient. 



Figure 4.4: Reconstruction results for noisy pressure data. Pressure data with 5% 
added noise (left; detector location varies in the horizontal and time in the vertical direction), 
reconstructed absorption coefficient using our single-stage approach (middle) the classical two- 
stage approach (right). 

Finally, in order to investigate the stability of the derived algorithms with respect to noise, 
we applied the single-stage and the two-stage algorithm after adding Gaussian white noise to the 
data with standard deviation equal to 5% of the maximal absolute data values. Note that for 
both, the single-stage and the two-stage algorithm, noise has only been added to the acoustic 
data. The reconstruction results for noisy data are shown in Figure 4.4. As can be seen both 
algorithms are quite stable with respect to data perturbations. However, again, the single-stage 
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approach yields better results and less artifacts than the two-stage algorithm. 


5 Conclusion 

In this paper we proposed a single-stage approach for quantitative PAT. For that purpose we derive 
algorithms that directly recover the optical parameters from the measured acoustical data. This is 
in contrast to the usual two-stage approach, where the absorbed energy distribution is estimated 
in a first step, and the optical parameters are reconstructed from the estimated energy distribution 
in a second step. Our single-stage algorithm is based on generalized Tikhonov regularization and 
minimization of the Tikhonov functional by the proximal gradient algorithm. In order to show 
that Tikhonov regularization is well-posed and convergent we analyzed the stationary radiative 
transfer equation (1.1), (1.2) in a functional analytic framework. For that purpose we used recent 
results of [ 20 ] that guarantees the well-posedness even in the case of voids. 

We presented results of our initial numerical studies using a simple limited angle scenario, 
where the scattering coefficient is assumed to be known. In this situation our single-stage algo¬ 
rithms has led to less artifacts than the two-stage procedure. More detailed numerical studies 
will be presented in future work. In that context, we will also investigate the use of multiple 
illuminations and multiple wavelength, which allows to also reconstruct the in general unknown 
scattering coefficient and Griineissen parameter. We further plan to investigate the use of more 
general regularization functionals such as the total variation in combination with the single-stage 
approach. 


Acknowledgements 

This work has been supported by the Tyrolean Science Fund (Tiroler Wissenschaftsfonds), project 
number 153722. We want to thank the referees for useful comments which helped to improve our 
manuscript, and for bringing reference [5( ] to our attention. 


References 

[1] M. Agranovsky, P. Kuchment, and L. Kunyansky. On reconstruction formulas and algorithms for the 
thermoacoustic tomography. In L. V. Wang, editor, Photoacoustic imaging and spectroscopy , chapter 8, 
pages 89-101. CRC Press, 2009. 

[2] H. Ammari, E. Bossy, V. Jugnon, and H. Kang. Reconstruction of the optical absorption coefficient 
of a small absorber from the absorbed energy density. SIAM J. Appl. Math., 71(3):676-693, 2011. 

[3] S. R. Arridge. Optical tomography in medical imaging. Inverse Prohl., 15(2):R41-R93, 1999. 

[4] H. Attouch, J. Bolte, and B. F. Svaiter. Convergence of descent methods for semi-algebraic and tame 
problems: proximal algorithms, forward- backward splitting, and regularized Gauss-Seidel methods. 
Math. Program., 137(1-2):91—129, 2013. 

[5] G. Bal, A. Jollivet, and V. Jugnon. Inverse transport theory of photoacoustics. Inverse Prohl., 
26:025011, 2010. 

[6] G. Bal and K. Ren. Multi-source quantitative photoacoustic tomography in a diffusive regime. Inverse 
Prohl, 27(7):075003, 20, 2011. 

[7] P. Beard. Biomedical photoacoustic imaging. Interface focus, 1(4):602-631, 2011. 


20 


[8] J. Bolte, S. Sabach, and M. Teboulle. Proximal alternating linearized minimization for nonconvex and 
nonsmooth problems. Math. Program., 146(1-2, Ser. A):459-494, 2014. 

[9] P. Burgholzer, J. Bauer-Marschallinger, H. Grim, M. Haltmeier, and G. Paltauf. Temporal back- 
projection algorithms for photoacoustic tomography with integrating line detectors. Inverse Probl., 
23(6):S65-S80, 2007. 

[10] J. Chen and Y. Yang. Quantitative photo-acoustic tomography with partial data. Inverse Probl., 
28(11): 115014, 2012. 

[11] E. Chouzenoux, J.-C. Pesquet, and A. Repetti. Variable metric forward-backward algorithm for 
minimizing the sum of a differentiable function and a convex function. J. Optim. Theory Appl., 
162(1):107-132, 2014. 

[12] P. L. Combettes and J.-C. Pesquet. Proximal splitting methods in signal processing. In Fixed-point 
algorithms for inverse problems in science and engineering, pages 185-212. Springer, 2011. 

[13] P. L. Combettes and V. R. Wajs. Signal recovery by proximal forward-backward splitting. Multiscale 
Model. Sim., 4(4):1168-1200, 2005. 

[14] B. Cox, J. G. Laufer, S. R. Arridge, and Paul C. Beard. Quantitative spectroscopic photoacoustic 
imaging: a review. J. Biomed. Opt., 17(6):0612021, 2012. 

[15] B. T. Cox, S. A. Arridge, and P. C. Beard. Gradient-based quantitative photoacoustic image recon¬ 
struction for molecular imaging. Proc. SPIE, 6437:64371T, 2007. 

[16] B. T. Cox, S. R. Arridge, P. Kostli, and P. C. Beard. Two-dimensional quantitative photoacoustic im¬ 
age reconstruction of absorption distributions in scattering media by use of a simple iterative method. 
Appl. Opt., 45(8):1866-1875, 2006. 

[17] B. T. Cox, T. Tarvainen, and S. R. Arridge. Multiple illumination quantitative photoacoustic to¬ 
mography using transport and diffusion models. In G. Bal, D. Finch, P. Kuchment, J. Scliotland, 
P. Stefanov, and G. Uhlmann, editors, Tomography and Inverse Transport Theory, volume 559 of 
Contemporary Mathematics, pages 1 12. AMS, 2011. 

[18] R. Dautray and J. Lions. Mathematical analysis and numerical methods for science and technology. 
Vol. 6. Springer-Verlag, Berlin, 1993. 

[19] A. De Cezaro and T. F. De Cezaro. Regularization approaches for quantitative photoacoustic tomog¬ 
raphy using the radiative transfer equation, http: //arXiv: 1307.3201, 2013. 

[20] H. Egger and M. Schlottbom. An L p theory for stationary radiative transfer. Appl. Anal., 93(6):1283- 
1296, 2014. 

[21] H. Egger and M. Schlottbom. Numerical methods for parameter identification in stationary radiative 
transfer. Comput. Optim. Appl., 2014. online first. 

[22] H. Egger and M. Schlottbom. Stationary radiative transfer with vanishing absorption. Math. Models 
Methods Appl. Sci., 24(5):973-990, 2014. 

[23] L. C. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics. Amer. 
Math. Soc., Providence, RI, 1998. 

[24] F. Filbir, S. Kunis, and R. Seyfried. Effective discretization of direct reconstruction schemes for 
photoacoustic imaging in spherical geometries. SIAM J. Numer. Anal., 52(6):2722—2742, 2014. 

[25] D. Finch, M. Haltmeier, and Rakesh. Inversion of spherical means and the wave equation in even 
dimensions. SIAM J. Appl. Math., 68(2):392—412, 2007. 

[26] D. Finch, S. K. Patch, and Rakesh. Determining a function from its mean values over a family of 
spheres. SIAM J. Math. Anal., 35(5):1213-1240, 2004. 


21 


[27] D. Finch and Rakesh. The spherical mean value operator with centers on a sphere. Inverse Probl., 
23(6):37-49, 2007. 

[28] J. Frikel and E. T. Quinto. Artifacts in incomplete data tomography - with applications to photoa¬ 
coustic tomography and sonar. arXiv: 1407.3453 [math.AP], 2014. 

[29] M. Haltmeier. Inversion of circular means and the wave equation on convex planar domains. Comput. 
Math. Appl, 65(7):1025-1036, 2013. 

[30] M. Haltmeier. Universal inversion formulas for recovering a function from spherical means. SIAM J. 
Math. Anal., 46(l):214-232, 2014. 

[31] M. Haltmeier, T. Schuster, and O. Sclierzer. Filtered backprojection for thermoacoustic computed 
tomography in spherical geometry. Math. Method. Appl. Sci., 28(16):1919-1937, 2005. 

[32] F. John. Partial Differential Equations, volume 1 of Applied Mathematical Sciences. Springer Verlag, 
New York, fourth edition, 1982. 

[33] G. Kanschat. Solution of radiative transfer problems with finite elements. In Numerical methods in 
multidimensional radiative transfer, pages 49-98. Springer, Berlin, 2009. 

[341 R. Kowar. On time reversal in photoacoustic tomography for tissue similar to water. SIAM J. Imaqinq 
Sci., 7(l):509-527, 2014. 

[35] R. A. Kruger, K. K. Kopecky, A. M. Aisen, Reinecke D. R., G. A. Kruger, and W. L. Kiser. Ther¬ 
moacoustic CT with Radio waves: A medical imaging paradigm. Radiology, 200(l):275-278, 1999. 

[36] R. A Kruger, P. Lui, Y. R. Fang, and R. C. Appledorn. Photoacoustic ultrasound (PAUS) recon¬ 
struction tomography. Med. Phys., 22(10):1605-1609, 1995. 

[37] P. Kuchment and L. A. Kunyansky. Mathematics of thermoacoustic and photoacoustic tomography. 
Eur. J. Appl. Math., 19:191 -224, 2008. 

[38] L. A. Kunyansky. Explicit inversion formulae for the spherical mean Radon transform. Inverse Probl., 
23(l):373-383, 2007. 

[39] A. V. Mamonov and K. Ren. Quantitative photoacoustic imaging in radiative transport regime. 
Comm. Math. Sci., 12(2):201 234, 2014. 

[40] M. Mokthar-Kharroubi. Mathematical topics in neutron transport theory. World Scientific, 1997. 

[41] W. Naetar and O. Sclierzer. Quantitative photoacoustic tomography with piecewise constant material 
parameters. SIAM J. Imaging Sci., 7(3):1755-1774, 2014. 

[42] V. P. Palamodov. Remarks on the general Funk-Radon transform and thermoacoustic tomography. 
Inverse Probl. Imaging, 4(4):693-702, 2010. 

[43] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. Experimental evaluation of reconstruction 
algorithms for limited view photoacoustic tomography with line detectors. Inverse Probl., 23(6):S81- 
S94, 2007. 

[44] K. Ren, H. Gao, and H. Zhao. A hybrid reconstruction method for quantitative PAT. SIAM J. 
Imaging Sci., 6(l):32-55, 2013. 

[45] A. Rosenthal, V. Ntziachristos, and D. Razansky. Acoustic inversion in optoacoustic tomography: A 
review. Curr. Med. Imaging Rev., 9(4):318-336, 2013. 

[46] A. Rosenthal, D. Razansky, and V. Ntziachristos. Fast semi-analytical model-based acoustic inversion 
for quantitative optoacoustic tomography. IEEE Trans. Med. Imag., 29(6):1275-1285, 2010. 

[47] T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge. A gradient-based method for quantitative 
photoacoustic tomography using the radiative transfer equation. Inverse Probl., 29(7):075006, 2013. 


22 



[48] O. Sclierzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen. Variational methods in 
imaging , volume 167 of Applied Mathematical Sciences. Springer, New York, 2009. 

[49] T. Schuster, B. Kaltenbacher, B. Hofmann, and K. S. Kazimierski. Regularization methods in Banach 
spaces, volume 10 of Radon Series on Computational and Applied Mathematics. Walter de Gruyter, 
Berlin, 2012. 

[50] N. Song, C. Deumie, and A. Da Silva. Considering sources and detectors distributions for quantitative 
photoacoustic tomography. Biomed. Opt. Express, 5(ll):3960-3974, 2014. 

[51] P. Stefanov and G. Uhlmann. Thermoacoustic tomography with variable sound speed. Inverse Probl., 
25(7):075011, 16, 2009. 

[52] T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. A. Arridge. Reconstructing absorption and scattering 
distributions in quantitative photoacoustic tomography. Inverse Probl., 28(8):084009, 2012. 

[53] K. Wang and M. Anastasio. Photoacoustic and thermoacoustic tomography: Image formation princi¬ 
ples. In Handbook of Mathematical Methods in Imaging, chapter 18, pages 781-815. Springer, 2011. 

[54] L. V. Wang. Multiscale photoacoustic microscopy and computed tomography. Nat. Photonics, 
3(9):503-509, 2009. 

[55] M. Xu and L. V. Wang. Photoacoustic imaging in biomedicine. Rev. Sci. Instrum., 77(4):041101, 
2006. 

[56] Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment. Reconstructions in limited-view thermoa¬ 
coustic tomography. Med. Phys., 31(4):724-733, 2004. 

[57] L. Yao, Y. Sun, and J. Huabei. Transport-based quantitative photoacoustic tomography: simulations 
and experiments. Phys. Med. Biol., 55(7):1917 1934, 2010. 

[58] R. J. Zemp. Quantitative photoacoustic tomography with multiple optical sources. Appl. Opt., 
49(18):3566-3572, 2010. 

[59] X. Zhang, W. Zhou, X. Zhang, and H. Gao. Forward-backward splitting method for quantitative 
photoacoustic tomography. Inverse Probl., 30(12):125012, 2014. 


23 



